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O ! ABSTRACT 

The Balloon-borne Large Aperture Submillimeter Telescope (BLAST) has made one square degree, 
deep, confusion limited maps at three different bands, centered on the Great Observatories Origins 
Deep Survey South field. By calculating the covariance of these maps with catalogs of 24 /im sources 
from the Far-Infrared Deep Extragalactic Legacy Survey (FIDEL), we have determined that the total 
I/-) '. submillimeter intensities are 8.60 ± 0.59, 4.93 ± 0.34, and 2.27 ± 0.20 nWm~ 2 sr" 1 at 250, 350, and 

£NJ ' 500 /Ltm, respectively. These numbers are more precise than previous estimates of the cosmic infrared 

background (CIB) and are consistent with 24 /xm-selected galaxies generating the full intensity of 
the CIB. We find that the fraction of the CIB that originates from sources at z > 1.2 increases 
with wavelength, with 60% from high redshift sources at 500 /im. At all BLAST wavelengths, the 
| relative intensity of high-z sources is higher for 24/im-faint sources than it is for 24/im-bright sources. 

Galaxies identified as active galactic nuclei (AGN) by their Infrared Array Camera (IRAC) colors are 
*-£h ■ 1.6-2.6 times brighter than the average population at 250-500 /im, consistent with what is found for 

»pHj X-ray-selected AGN. BziGselected galaxies are found to be moderately brighter than typical 24 fim- 

q selected galaxies in the BLAST bands. These data provide high precision constraints for models of 

the evolution of the number density and intensity of star forming galaxies at high redshift. 

Subject headings: cosmology: observations — cosmology: diffuse radiation — submillimeter — galax- 
ies: evolution — galaxies: starburst 
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1. INTRODUCTION 



A decade ago, a uniform cosmic infrared background 
radiation (CIB) was discovered in Far Infrared Abso- 
lute Spectrophotometer (FIRAS) d ata from the Cosmi c 
Background Explorer ( COBE) by iPuget et all (|1996f) . 
and later confirmed by iFixsen et alT i 19981 ). The back- 
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ground, which peaks in intensity at a wavelength of 
w 200 um, is as bright as all the optical light in the Uni- 
verse (Hauser & Dwck 2001). It is presumed to be ther- 
mal re-radiation of optical and UV starlight absorbed 
by dust grains. At around the same time as the dis- 
covery of the CIB, observations with the Submillimetre 
Com mon User Bolometer Array (SCUBA, Hollan d" et al.l 
1999) at the James Clerk Maxwell Telescope revealed a 
population of dusty high redshift starburst ga laxies form- 
ing stars at a rate of ~ 1000 per year (jBlain et alJ 
[TflM [Hughes et alJ[l99l IBareer et al.lll998l K 

Over the ensuing decade, a substantial effort has been 
devoted to understanding the implications of this back- 
ground radiation, and the task has not been easy. Only 
the few very brightest submillimeter galaxies are vis- 
ible with SCUBA in blank-fields above the confusion 
limit, and these compr ise a small fraction of the to- 
tal CIB at 850/xm (e.g lHughes et al l 119981 : iScott et all 
l2002HBorvs et al.ll2003t iCoppin et al.ll2006l) . Full under- 
standing relics on statistical measures of submillimeter 
galaxies too weak to be detected individually. Attempts 
to probe the contribution to the CIB by faint sources 
have been made by looking at objects lensed by clus- 
ters flCpwie et a l. 2002; Sm ail et al.ll2002l : iKnudsen et alJ 
I2008f) . iKnudsen et all (|2008f) have resolved 100% of the 
850 /im background by reaching unlensed flux density 
limits of 0.1 mjy, but the reliability of the lensing inver- 
sion introduces uncertainty in the estimate, and direct 
detection of these sources is desirable. Furthermore, the 
CIB signal peaks at wavelengths where the atmosphere 
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is opaque and has dropped by almost two orders of mag- 
nitude before it reaches the reliable atmospheric window 
SCUBA exploits. The bulk of the CIB is not visible from 
the ground. 

Since both the numerical density and the emission 
spectra of galaxies evolve with cosmic epoch, the CIB 
is a convolution of these evolution functions for a wide 
range of galaxy types. Because of this, the shape of the 
CIB contains a wealth of information about the history of 
star formation. However, working from measurements at 
only one wavelength, the convolution can not be inverted 
to extract an evolution function. One expects a different 
mixture of galaxy types, and a different redshift distribu- 
tion, to emerge from observations of the CIB at different 
wavelengths. Measurements at several wavelengths are 
therefore required to constrain models of these underly- 
ing distributions. 

BLAST, the Balloon-borne Large Aperture Submil- 
limeter Telescope, was designed to carry out the program 
of characterizing the CIB over a range of wavelengths 
near to its peak. We report here on the success of that 
program. BLAST has produced deep, confusion-limited 
maps at three wavelengths (250, 350, and 500 /an), where 
the CIB produces substantial flux density, and we have 
used a uniform and c arefully constructed catalog of 
24 /zm-selected galaxies (jMagnelli et al.l [20091) whose in- 
tensities we can find statistically in the BLAST data. 
These data bridge the gap between similar analyses at 
shorter wavelengths, using data from the Spitzer Space 



i gtns, using data t 
Telescope (e.g. IDole et alJ l2006f L an d at longer wave- 
lengths, using data from SCUBA (e.g. IWang et alj|2006t 
iSerieant et~aLl|2008fl. 

iDevlin et al.l (|2009T ) established that the full intensity 
of the CIB is resolved statistically into flux density pro- 
duced by identifiable 24 /im-selected galaxies. In this pa- 
per we examine that relation in more detail by dividing 
the 24 /zm c atalogs by brightness and color. A compan- 
ion paper bv lPascale et al l (|2009l ) uses spectroscopic and 
photometric redshifts to constrain the star formation rate 
history. Chapin et al. (in prep.) will explore the implica- 
tions of these new BLAST results on our understanding 
of models of galaxy evolution. 

We analyze these high-S/N BLAST maps of confusion 
using techniques that differ from those appropriate for 
catalogs of isolated point sources; we use the covariances 
between our maps and external catalogs, a technique 
know n as "stacking " , to m easu re the background . Sim- 
ilarly, IDevlin et all (|2009f ) and iPatanchon et all (|2009l ) 
use a a P(D)" fluctuation analysis, rather than counts 
of individual sources, to provide the most accurate 
measurements of the underlying source counts, while 
iViero et al.l (|2009l ) study correlations in the background. 
The possibility of determining number counts, spectral 
energy distributions, and clustering biases directly from 
correlations in the maps without the requirement to first 
extract a c a talog of point sources was pointed out by 
iKnox et all (|200l1 ). 

The layout of this paper is as follows: in §[2]thc BLAST 
observations and external data are described. In §[3] we 
develop the stacking formulae for measuring the contri- 
bution to the total surface brightness detected in a map 
produced at positions from an external catalog. Finally, 
in §|4] we use 24 (im and optical catalogs, and subsets 
based on color cuts, to determine the total CIB mea- 



sured by BLAST, as well as the relative contributions of 
low- and high-redshift sources, AGN, and BzK galaxies. 
Additionally, we include an analysis of deep extragalac- 
tic SCUBA 85 jitm o bservations, similar to those used 
bv IWang et all ([200l . 

2. OBSERVATIONS 
2.1. BLAST 

BLAST is a 1.8 m stratospheric balloon-borne tele- 
scope that operates at an altitude of approximately 
35 km, above most of the opaque atmospheric water va- 
por that renders observations at all but a few narrow 
submillimeter bands from the ground impossible. Obser- 
vations are undertaken simultaneously with three broad- 
band bolometric imaging arrays with central wavelengths 
250, 350, and 500 /im; this camera is a prototype of the 
Spectral a nd Photometric Im aging Receiver (SPIRE) for 
Herschel (IGriffin et al.ll200l . 

We use data from the 11-day BLAST flight in 
2006 from McMurdo Station, Antarctica. The undcr- 
illuminatcd BLAST primary produced nearly diffraction- 
limited beams with full-width at half-maxima (FWHM) 
of 36", 42", and 60" in each band, respectively. This suc- 
cessful flight produced significantly deeper- and wider- 
area extra-galactic maps than existing 350 and 450 /im 
ground-based observations, and the first 250 /im maps, 
near the peak in the CIB. Deep and wide blank- field 
surveys were conducted in the Great Observatories Ori- 
gins Deep Survey South (GOODS-S) field, which we la- 
bel BLAST GOODS-S Deep (BGS-Deep) and BLAST 
GOODS-S Wide (BGS-Wide), respectively. These maps 
are centered on the Ch andra Deep Field So uth. A cov- 
erage map is shown in Pascalc et alJ (|2009l ) . A second 
intermediate-depth field near to the South Ecliptic Pole, 
which will be the subject of a separate study, was also 
observed. In addition, several low-redshift clusters and 
high-redshift radio galaxies were targeted to sample bi- 
ased star-forming regions of the Uni verse. Further detail s 
on the instrument may be found in lPascale et al.l ((2008) , 
and the flight perfor mance and calibrat ion for the 2006 
flight are provided in iTruch et alJ (|2009l ). 

In this paper we focus on the deepest BLAST maps 
of BGS-Deep 15 wh ich cover an area of approximately 
0.8 square degrees (jDevlin et al.l 120091 ) . and completely 
encompass the Extended Chandra Deep Field South 
(ECDF-S), which in turn encompasses the smaller 
GOODS-S and Hubble Ultra Deep Field South fields. 
The maps in all BLAST bands are confusion limited, such 
that instrumental noise itself contributes only a fraction 
(~ 50%, see below) to the r.m.s. of the map, and therefore 
provide high S /N measurements of the spatial variations 
in the intensity of the CIB. 

All BLAST time-stream detector data are reduced us- 
ing a common pipeline to identify spikes, correct time- 
varying detector responsivities, and deconvolve the lag 
induced by thermal time constants for the bolometers. 
Maps are produced from a combination of these cleaned 
data with the pointi ng solution using a maximum- 
likelihood algorithm (jPatanchon et all 120081 ). The ab- 
solute calibration is based on regular observations of the 
evolved star VY CMa, which results in uncertainties of 

15 The BLAST maps used in this paper are available for down- 
load at http://blastexperiment.info/ 
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Pixel signal — to-noise ratio 

Fig. 1. — The distribution of the pixel S/N ratio in BGS-Deep 
map at 500 (histogram). The scale on the y-axis is the number 
of pixels per bin of width 0.1. The positive tail is due to bright 
sources in the map. Because the map has a mean of zero, the peak 
of the distribution is shifted to the left of zero, as indicated by 
the vertical dotted line. A Gaussian (thick solid curve) is fit to 
the negative side of the histogram. For reference, the measured 
instrument noise (a Gaussian, with width 1.0 in units of S/N) is 
over-plotted (thin dashed line). The map is clearly dominated by 
confusion noise. 

approx imately 10% that a re strongly correlated between 
bands (|Truch et all 12009ft . The resulting maps may be 
thought of as the optimal weighting of the data on all 
spatial scales such that both point sources and diffuse 
structures are accurately reproduced within the limita- 
tions of the data. However, while all spatial frequencies 
(inverse spatial scales) from the mean level up to the 
Nyquist frequency are estimated, the largest scales (in- 
cluding the absolute brightness) are not well constrained 
due to detector drift and other systematics. We have 
therefore suppressed all scales larger than approximately 
9', 7.'5, and 8' at 250, 350, and 500 /mi, respectively - 
this procedure explicitly sets the mean of each map to 
zero. 16 

We show in Figure [T] the noise properties of the 500 /im 
BGS-Deep map. We plot the distribution of the ratio of 
map pixel values to instrumental noise, cxinst, estimated 
by the map-maker, which propagates uncertainties in the 
time-domain detector noise power spectra. The distribu- 
tion is well-described by a Gaussian with an excess at 
positive flux density, due to bright sources in the map. 
The shape of the histogram is due to both instrumen- 
tal noise and confusion, cr con f. We estimate the confu- 
sion noise by subtracting the instrumental noise from the 
map r.m.s. in quadrature, a con { = (c^ap — ^inst) 1 ^ 2 - The 
map r.m.s. is 16mJy and the median instrument noise is 
6.6 mjy, thus we conclude that the noise due to source 
confusion is CT C onf = 15mJy at 500 /im. Similar analy- 
ses give confusion noises of 17 and 21mJy based on in- 
strumental noises of 8.6 and llmJy at 350 and 250 /im, 
respectively. The maps are clearly dominated by con- 
fusion; this consists of a high S/N measurement of the 
shape produced by galaxies much too close together to 
be resolved individually. 

16 The filter is, in practice, anisotropic, with greater suppres- 
sion of the more poorly constrained modes orthogonal to the scan 
direction. 



2.2. External Catalogs and Data 

In order to estimate the contribution of sources in the 
BLAST maps to the CIB, we stack the map at positions 
from catalogs measured at other wavelengths. In the 
following subsections we describe these external catalogs. 
The details of the stacking method are described in § 13 11 

2.2.1. FIDEL /SIMPLE 

Recently it has been shown that deep Multiband Imag- 
ing Photometer for Spitzer (MIPS) 24 /tm catalogs con- 
tain significant fractions of the sources tha t produce the 
CIB at 70 and 1 60 /tm (|Dole et all 120061 ) and 850 /tm 
(jWang et al.ll2006P ). wavelengths that bracket the BLAST 
coverage. The deep est 24 /tm cata l og in BGS-Deep 
has been produced by[Magnelli et alJ ()2009| ). combining 
GOODS-S with Far-Infrared Deep Extragalactic Legacy 
survey (FIDEL) data that cover a total of 0.21 square 
degrees. The FIDEL catalog was constructed using a 
higher-resolution Infrared Array Camera (IRAC) cata- 
log as a positional prior. This deep catalog from SIM- 
P LE (Spitzer IRAC / MUSYC Public Legacy in ECDF- 
S, Gawiser et al. 2006) enabled the deblending of sources 
with separations as small as 0.5FWHM in the MIPS 
map, resulting in an extremely faint flux density limit 
5*24 15 /tJy. The faintest source in the catalog has 
.SW = 13 /tJy, but the detections at 5*24 < 30 /tJy should 
be considered tentative. The effects of this are discussed 
in S l4~5l 

Furthermore, since the catalog is a 24 /tm-detected sub- 
set of SIMPLE, there also exist IRAC 3.6, 4.5, 5.8, and 
8.0 /im flux densities for each object. The depth, relia- 
bility, and wavelength coverage of this catalog have been 
essential to this work. 

2.2.2. MUSIC 

Unlike the FIDEL catalog which has a fairly simple 
mid-IR selection function, the MUltiwavelength South- 
ern Infrared Catalog (MUSIC) uses a near-IR selection 
with a much more heterogeneous data set combining 
ACS and IRAC maps with ground-based JHK Very 
Large Telescope (V LT) data and optical spectroscopy 
(jGrazian et al.ll200l . The primary use of this catalog 
is to obtain spectroscopic and optical photometric red- 
shifts, although we also check the contribution to the 
CIB using stacking. 

2.2.3. SCUBA 

We extend the measurements made here with BLAST 
to lower frequencies using SCUBA data. The largest 
deep extragalactic SCUBA observations are in the 
GOODS North (GOODS-N) field. Data from this field 
have already been used to measure the 850 /tm c ontri- 
bution to the CIB (|Wang et all [20061 IPopl [2001 . but, 
for consistency, we perform our own analysis using tech- 
niques identical to what we have done with the BLAST 
data. The map used here was produced by iBorvs et alJ 
(2003), combin ing data from se veral different groups, and 
extended by iPope et all (|2005f ). covering an area of ap- 
proximately 0.06 square degrees. 17 The noise in these 
maps is not uniform, varying from approximately 0.5 

17 The SCUBA map o f GOODS-N us ed in this paper i s available 
at http://www.noao.edu/staff/pope/DATA/SCUBA.html 
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to 4mJy. The existin g Spitzer catalogs in GOODS-N 
(|Papovich et al.l 120041 ) are of similar depths to the FI- 
DEL catalogs of GOODS-S. The size of the SCUBA and 
24 fim overlap region in GOODS-N is 0.055 square de- 
grees. 

3. METHODOLOGY 
3.1. Stacking 

Determining the mean flux density of a population of 
sources that are individually too dim to be detected by 
examining their effect on a confusion limited map is not 
new. The main approach goes by the name of "stacking" . 
Because this work makes such heavy reliance on the tech- 
nique, and because technical questions arise about the 
generalization to very high source density or the wisdom 
of excluding bright sources, etc., we review the basis of 
the technique here. We find that many of these miscon- 
ceptions are avoided when one realizes that the technique 
is really one of taking the covariance of the map with the 
catalog. 

Imagine we have a map of the sky where D 3 is the 
flux density in each pixel j after convolution with the 
instrumental point-spread function. Suppose also that 
we have one or several independent catalogs of sources 
made from other experiments; catalog C a has N J a sources 
in pixel j, and we would like to know the mean flux 
density, S a , of the sources in C a . We denote the mean 
of Nl the average number of sources per pixel in 

list C a ■ If all sources in a catalog produce identical flux 
density, S a , then, along with whatever else is in the sky, 
there will on average be a contribution of S J a — S a N^ to 
each pixel. 

The BLAST maps have zero mean. We can write the 
flux density in the map as 

Dj = S 1 (N( - m) + S 2 (N} - M2 ) + ... + Wj , (1) 

where Wj is the contribution of detector noise in pixel j, 
and, strictly speaking, the S a form the complete set of all 
objects in the Universe. We have suppressed the mean by 
subtracting the S a /i a for each catalog from every pixel. 
We additionally require that Wj has a mean of zero. We 
imagine that the sources in the catalog are not correlated, 
such that N J a is a random, Poisson-distributed number. 
Furthermore, assume that no two lists are correlated, so 
that 

((Ni-fi a )(N^-^))=0, (2) 

Our goal is to determine the mean BLAST flux density 
per source in list C a from knowledge of the BLAST maps, 
Dj, and the locations, N^, of the sources in C a , but 
without any other information, and in particular without 
requiring any knowledge of other lists Cp (a ^ (3). 

Nj is a function that has a shape on the sky, and mea- 
suring the amplitude of this shape in the map determines 
S a , the mean source brightness. The covariance of Dj 
with Nj provides the maximum likelihood estimate 18 of 
this amplitude: 

Cov(D j ,Ni) = -}-^2D j Ni 

v 3 

18 This is the maximum likelihood estimate only if the map pixel 
noises Wj are uniform. This is not the case for the BLAST maps. 
We address this issue later on. 



= S'cCT^j , (3) 

where A p j x is the total number of pixels in the map, 
and we have dropped terms in N^Ni and N-^Wj, which 
vanish in the sum. The first equality in Equation holds 
because Dj has a mean of zero, and the last equality 
follows from the definition of variance. The variance of 
a Poisson distribution is 

so the net result is that the zero-lag cross-correlation 
(covariance) of a catalog with the map divided by the 
mean number of sources per pixel is an estimate of the 
average flux density per source, 

Sa = <Mg**a ( 5 ) 

A final re-arrangement facilitates the use of Equation[5j 
Notice that the sum in Equation [3] runs over all pixels, 
with the weight of each pixel proportional to the number 
of catalog sources found within it, and that zero weight 
is given to pixels that do not contain a catalog source 
(Nj = 0). This can be re- written as a sum, instead, over 
all catalog entries with unity weight, 

1 v pix /^a ' fc a , 

3 & 

where k is the index of sources in catalog C a , Dk is the 
measured flux density in the map pixel that contains the 
fc th catalog entry, and n a is the total number of cata- 
log entries, n a — N p i x ^i a . This extraordinarily simple 
formula can be used to probe the properties of sources 
much too crowded to be detected individually, and also 
those with flux densities that are much fainter than the 
typical thresholds of source catalogs derived only from 
the map itself. Since the map pixel noises Wj are not 
uniform across the map, we weight the mean in Equa- 
tion [5] by the inverse pixel variance to maximize the S/N 

ratio of S a . We explicitly subtract the weighted means 
of the BLAST maps. 

Equation [6] provides a robust estimate of the mean 
brightness per source even when there are other, possibly 
substantial, contributors to the flux density present, Cp. 
This is provided that iV a is Poisson distributed, and N£ 
is not correlated with either the detector noise or sources 
in Cp. In other words, the effect of other sources on the 

estimator S a is to provide an additional source of noise. 
This noise may potentially be asymmetric, but it has a 
mean of zero, such that S a is unbiased. 

Similarly, a catalog C a can be subdivided into disjoint 
subsets, and the mean brightness due to each subset can 
be measured without bias. We use this fact to split up 
our catalogs based on IRAC colors and brightness at 
24 fim. 

We are now in a position to address the proper han- 
dling of sources that are bright enough to be easily recog- 
nized in the maps, for example the sources in a BLAST 
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5 o~ catalog. We have shown that S a , our estimate of S a , 
is not affected by either the presence or the removal of 
flux density from other source lists C73 that are uncor- 
rected with C a . However, since the sum of confusion 
noise and detector noise, S a N^ + Wj, will cause sources 
near the threshold to be accidentally included or ex- 
cluded from the BLAST catalog, any list made from the 
BLAST maps themselves will be artificially correlated 
with all the terms in Equation [TJ Furthermore, since the 
BLAST-generated bright source catalog depends on the 
sum of the other terms in Equation [TJ excision of the 
flux density from such a catalog will artificially correlate 
the remaining terms, such as (JVjJ — fJ, a ) and Wj. This 

introduces a bias in our estimator S a that is difficult to 
quantify. In all the following work, stacking is performed 
on the full BLAST maps, including any bright sources 
they contain. 

We reiterate that that this formulation for stacking is 
slightly di fferent than that u sed by other authors. In 
particular, IDole et aTl (|2006T ) performed aperture pho- 
tometry on their stacked maps, effectively subtracting 
a local background. T he stacking analysis performed 
bv lPascale et al.l (120091) uses a similar technique to that 
of IDole et al l (|2006f ). and we have confirmed that their 
stacking calculations give results consistent with those 
presented here. 19 

The relation between the average BLAST flux density 
per source, S„, and the background specific intensity, /„, 
is S^n/fl, where v is the BLAST frequency and n/fl is 
the number of sources per solid angle. Note that the in- 
tensity of the background is often defined through the 
product vl v , which is then related to vS v . In what fol- 
lows, we use / = vl v . 

Near to its faintest limit, any catalog is only partially 
complete, so we estimate n/fl as n c /(Cf2) where n c is the 
number of catalog entries, C is the completeness of the 
catalog, assumed to be a function only of flux density, 
and O is the solid angle of the catalog. Completeness 
of the FIDEL 24 catalog is measured by comparing 
the differential flux density distribution of the catalog to 
publis hed completeness co rrected counts. Above 400 ^Jy 
we use IShupe et all (120081 Table 2) and b etween 35 p Jy 
and 400 /zJy we use lPapovich et all (|2004 Table 2) after 
correcting their flux densities for a 5.6% difference in cal- 
ibration that arises from choices in Spitzer photometry 

methods. The smooth function C = [l + (A/S24) 13 ] , 
where A = 16/xJy and (3 — 1.8, interpolates this correc- 
tion. This correction is C = 0.56 at our lowest 24 /im 
flux density bin and has a small effect on our final result 
(see S ET~2"1) . 

3.2. Uncertainties 

To estimate the uncertainty of Equation[S]algebraically 
for a catalog C a , one would need to know the scatter pro- 
duced by its complement — the catalog of all sources not 
in C a that contribute to the background (in addition to 
sources of instrumental noise). In practice, the comple- 
ment is not known, so we establish the uncertainties and 
possible biases of our measurements by generating ran- 
dom catalogs and stacking them on the actual BLAST 

19 Pascalc ct al. (2009) use a slightly different reduction and 
spatial selection of the BLAST data than we use here, so we don't 
expect the results to be exactly the same. 



maps under analysis. We find, as expected, that the un- 
certainties are Gaussian-distributed and scale as the map 
r.m.s. (including confusion noise) divided by the square 
root of the number of catalog entries. See Figured! and 
notice the high precision of achieving a null, i.e. zero in- 
ferred average flux density, from stacking on a catalog 
of random positions. The uncertainties arise from sam- 
pling the maps at random locations; since the confused 
flux density in the three separate BLAST maps is cor- 
related, the uncertainties for stacking are correlated be- 
tween BLAST bands. We measure the correlations from 
the simulations and find P12 = 0.70, P13 = 0.68, and 
P23 = 0.70, where p a b is the Pearson correlation coeffi- 
cient between bands a and b, and 1, 2, and 3 refer to 250, 
350, and 500 /im, respectively. 

3.3. Limitations of the method 

The stacked signals only return the mean source flux 
density if the catalog is Poisson distributed, since we have 
used the fact that the variance of the number of catalog 
entries in a region of a given size equals the mean. When 
this assumption is violated, the relation between the co- 
variance of a map with a catalog and the corresponding 
mean flux density becomes complicated. If the catalog 
C a is clustered at the scale of the BLAST beams, it is 
very easy to overestimate the mean flux density S a and 
therefore overestimate the contribution of a given catalog 
to the total sky intensity. 

This problem is strikingly apparent when we per- 
form the covariance with the MUSIC catalog, which has 
~ 18,000 sources in 0.035 square degrees (500,000 sources 
per square degree). We find that the background inten- 
sity inferred by stacking the BLAST maps on the MUSIC 
catalog exceeds the FIRAS values by factors of 2.5-3.0. 
We suggest that clustering in the catalog (whether due 
to selection biases or real clustering on the sky) is the 
cause of this overestimate. 

The level of clustering that is measurable in a given 
area is strongly dependent on the source density, since, 
for a given level of clustering, the accuracy with which 
the mean can be measured goes as the square root of the 
number of sources in that area. We show here that the 
MUSIC catalog is strongly clustered at the scale of the 
BLAST beams and thus we should expect a biased mea- 
surement of the background, but that the FIDEL catalog, 
with ~9100 sources in 0.21 square degrees (43,000 per 
square degree), is well-behaved at the appropriate scales. 

We have measured how close to Poisson-distributed the 
MUSIC and FIDEL catalogs are. We place 500 circles of 
diameter D at random within the catalog area and count 
the number of catalog entries in each. Figure [3] shows the 
ratio of the variance to the mean as a function of circle 
diameter, for MUSIC and FIDEL. If the catalogs were 
uniform random distributions, this ratio would be unity. 
The MUSIC catalog shows a substantial excess variance 
at all angular scales above a few tens of arc seconds. It 
is not a surprise that galaxy locations are correlated and 
that MUSIC is deep enough to measure that. We also 
note that MUSIC is an extremely heterogeneous cata- 
log, consisting of a variety of pointed observations. In 
contrast, FIDEL is essentially a flux-limited catalog pro- 
duced from a nearly uniformly sampled map. We con- 
clude that on the angular scale of the BLAST beam sizes, 
one should anticipate that covariance with MUSIC will 
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Fig. 2. — Quantification of errors in FIDEL stacking measurements. We calculate the average BLAST flux density at N random positions 
in the each map, where N is equal to the number of sources in the FIDEL catalog. Upper row: histograms of 10 stacking measurements 
at random positions. The scale on the y-axis is the number of simulations per 5 (j,Jy bin. Lower row: difference between the histograms 
and Gaussians centered at 0.0 with width equal to the r.m.s. of the distributions — clearly the histograms are very well described by 
Gaussians. We use the r.m.s. of these distributions as the errors in the stacked background values. We note that the errors between bands 
are strongly correlated, with Pearson correlation coefficients p~0.7 (see Table[TJ in all cases. 

of many arc minutes. We have tested whether this differ- 
ence arises because FIDEL sources are intrinsically dif- 
ferent from MUSIC sources, or is instead simply a fea- 
ture of the shallower depth. We repeat the tests shown 
in Figure [3] for just the 1200 FIDEL sources that lie in 
the MUSIC region and for several random subsamples 
of 1200 MUSIC catalog entries. None of these curves 
is statistically distinguishable from unity at any angular 
scale, and we conclude that the Poisson variance associ- 
ated with sampling MUSIC at the FIDEL number den- 
sity dominates over the correlations in galaxy locations 
detected by MUSIC. 

4. DISSECTING THE SUBMILLIMETER BACKGROUND 

4.1. Total Intensity 

We stack the BLAST maps on the FIDEL 24 fim cata- 
log, using the methods described in § 13.11 Completeness- 
corrected results are shown in Figure 2] and listed in 
Table [D along wi th v alues of th e CIB measured by 
iFixsen et all (|1998f l and lDole etail (|2006| ). The BLAST 
measurements are consistent with having resolved 100% 
of the background, as measured by FIRAS. These num- 
bers should be seen as lower limits, though; see § 14.21 
Even including calibration uncertainty, the BLAST to- 
tal intensity values are twice as precise as FIRAS. In the 
rest of this paper, total intensity of the CIB refers to the 
BLAST values listed in column 2 of Table □ We note 
that while the completeness of the FIDEL catalog is un- 
certain (§ 13. 1|) , the correction to the BLAST background 
intensities is small (< 10%, see § 14. 2|) . 

We examine the mean BLAST flux density per source 
as a function of 24 /im flux density by dividing the FIDEL 




Diameter farcminl 



Fig. 3. — Clustering in the MUSIC and FIDEL catalogs as a 
function of scale. The ratio of the variance of the number of sources 
to the mean number of sources within a circle is plotted against the 
diameter of the circle for both the MUSIC and FIDEL catalogs. For 
a Poisson-distributed catalog, this ratio would be unity at all scales. 
Each point in the graph is the variance (divided by the mean) of the 
number of entries that fall within 500 randomly-placed circles of a 
given size within the respective survey area. The broken vertical 
lines indicate the BLAST FHWMs. Notice that the FIDEL catalog 
retains a Poisson distribution to larger scales than the BLAST 
beams while the much deeper and heterogeneous MUSIC catalog 
shows excess variance arising from the clustering of galaxies at all 
angular scales. 

provide a biased estimator of total intensity. 

The FIDEL catalog, which has a substantially lower 
surface density, does not show excess variance until scales 
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TABLE 1 

Corrected stacked intensities 







RT, AST 






Band 


Total 


Low-z 


High-z 


FIRAS 


Qua) 


(nWm- 2 sr- 1 ) 


(nWm- 2 sr- 1 ) 


(nWm^sr" 1 ) 


(nWm- 2 sr" 1 ) 


250 


8.60 ±0.59 (1.04) 


5.18 ±0.45 (0.69) 


3.42 ±0.37 (0.50) 


10.4 ±2.3 


350 


4.93 ±0.34 (0.68) 


2.44 ±0.26 (0.39) 


2.49 ±0.21 (0.37) 


5.4 ± 1.6 


500 


2.27 ±0.20 (0.36) 


0.89 ±0.15 (0.20) 


1.38 ±0.13 (0.22) 


2.4 ±0.6 



Note. — The quoted errors are measurement uncertainties only. The num- 
bers in parentheses are errors including calibration uncertainties. The errors be- 
tween the BLAST bands are strongly correlated because BLAST observes similar 
confusion-limited structure at all wavelengths. The Pearson correlation coefficients 
are (pi2, P13, P23) = (0.70,0.68,0.70) for measurement-only uncertainties, where 1, 2, 
and 3 refer to 250, 350, and 500 (im, respectively. The coefficients for the full uncer- 
tainties are (0.90, 0.88, 0.91), (0.86, 0.80, 0.82), and (0.84, 0.83, 0.90) for the total, low-z, 
and hig h-z: stacks, respectiv ely. We note that these numbers are slightly different from 
those in Devlin et al. (2009), due to a small update in the calculation methods. 
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Fig. 4. — The total intensity of the far-IR, background (CIB) 
in the BLAST bands associated with 24 fim galaxies is plotted 
against wavelength. The BLAST error bars include calibration 
uncertainties and are highly correlated between bands. The point 
at 850 pm labele d SCUBA is a simi lar st ack in GOODS-N using 
data described in lBorvs et~aTl H2003fl and lPope eTaLI (120051) . The 
hatched region marked "Fixsen" is the absolute spe ctrum of the 
CIB determined from FIRAS observations on COBE ( Fixsen ct al. 
ll99STy The hatched re gion at A < 160 urn labeled "Dole" is the 
background reported by Dole et al. (2006) using a similar stacking 
analysis of Spitzer data in GOODS-N. 



catalog into bins of S24 (Figured]). The brightest sources 
in the BLAST catalog are a factor of 10 times brighter 
than the mean brightnesses of 24 //m-selected sources 
shown in this figure; it is tempting to imagine that these 
curves extend to the right, but no brighter 24 fim sources 
exist. We conclude that the BLAST sources have anoma- 
lous ^rt.art 1 yet the full Sblast > 3<t catalogs 
pevlin et al.l 12009) comprise only 10-15% of the total 
intensity determined from stacking the full FIDEL cata- 
log, since the number density of these objects is low. 

The data do not allow a linear relation between flux 
density at 24 /im and flux density in the BLAST bands, 
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500 
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Fig. 5. — The mean BLAST flux density of FIDEL sources is 
plotted against 24 (im flux density for each of the three BLAST 
bands. The 350 and 500 (im curves have been displaced horizon- 
tally by —5 and ±5%, respectively, for visual clarity. The line near 

1/2 

the uppermost curve is Sblast <x S 24z , shown for reference. The 
vertical lines indicate the MIPS flux densities for which mean spec- 
tra are plotted in Figure [6] The number of FIDEL 24 /an sources 
per bin decreases from ~ 2600 at the faint end to 3 at the bright 
end. The lowest flux density bins are based on tentative 24 pm 
detections, and are therefore biased low (see § 14,51 1, 



such as any model of the form 

Ablast 



S 



BLAST 



A 



24 /im 



S, 



24' 



(7) 



with (3 = 1, would imply for any value of p. This sug- 
gests that the BLAST-detected light comes from a dis- 
tribution of sources with different spectral energy distri- 
butions (SEDs), or with a range of redshifts that varies 
significantly as a function of S24. This conclusion relies 
on the very large dynamic range, a factor of 500 in S24, 
available from the FIDEL catalog. Even for a power law 
relation, allowing (3 to differ from 1, the data do not sup- 
port a constant value of (3 except perhaps at high flux 
densities. Figure [5] shows how the 250 /mi curve diverges 
from (3 = 0.5 at low flux densities. 

Spectra of the mean flux density in BLAST bands of 
sources at three fixed 24 /im flux densities are shown in 
Figure[6] This plot shows two features very clearly: first, 
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Fig. 6. — Spectra of the mean flux density in BLAST bands of 
sources of fixed 24 fim flux density. The three curves correspond 
to the 24 fim flux densities indicated by the three vertical lines in 
Figure [5] The point at 24 (im on each curve is the central flux den- 
sity of the bin defining each sub-sample and the error bar indicates 
the width of that flux density bin. The circles at 70 and 160 fim 
on the lower two curves are the mean flux densities in those bins 
determined by a similar stacking analysis in Dole ct al. (2006, Fig- 
ure 7). 70, 160, and 850 /im flux densities are not available for the 
upper curve, so we plot dashed lines following the same slopes as 
the middle curve. A clear trend of increasing mean flux density and 
apparent temperature (spectrum peaking at shorter wavelengths) 
is observed in the BLAST bands with increasing 24 fim flux density. 

that the average BLAST flux density is positively corre- 
lated with 24 fim flux density; and second, that fainter 
24 fim sources appear cooler (as suggested by the shal- 
lower slopes through the BLAST wavelengths). This lat- 
ter point is probably due to the fact that higher-redshift 
sources have predominantly fainter 24 fim flux densities 
(see, for example, the strongly no n-Euclidean region o f 
the source counts at S24 < 1 mJy in lPapovich et aU2004D 
- if the average rest-frame galaxy dust temperature 
does not evolve appreciably, then the mean spectrum of 
fainter/more distant objects will undergo greater redshift 
and hence appear cooler. We note, however, that these 
curves are not well described by a single-object SED, and 
are certainly averages over a wide range of source types. 
The effects of a rough cut on redshift are presented in the 
next section, and detailed BLAST stacks as a function 
of redshift are explored in iPascale et al.l (|2009f) . 

4.2. Division in Redshift 

We use IRAC 3.6-8.0 ^m flux densities of the FIDEL 
sources to make a cut in a color-color plane to broadly 
classify the sources as either high- or low-redshift (Fig- 
ure [7|). The sources (small grey dots in the figure) mostly 
lie in two clouds, one in the lower-left quadrant and the 
other in the upper-left quadrant of the figure. For il- 
lustration, we over-plot the observed colors of a galaxy 
template at a range of lookback times (iit>) ; ranging from 
0-12 Gyr (black circles). The template is the SBc galaxy 
VCC 1987 from lDevriendt et all (|1999f ). chosen because 
it is publicly available and exhibits the major features of 
the color-color space occupied by the FIDEL catalog. We 
note that the clouds correspond to regions of color-color 
space where this template galaxy lingers for several Gyr. 
The track traced out by the local star-forming galaxy 
M82, thought to be representative of high-redshift sub- 
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-10 12 
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Fig. 7. — IRAC color-color plot for sources in the FIDEL cata- 
log (small grey dots). The sources lie in two partially-overlapping 
clouds. The black circles rep resent the colors of an SBc galaxy tem- 
plate (Dcvricndt ct al. 1999, VCC 1987) observed at a range of red- 
shifts. The circles are linearly spaced in lookback time; the small 
circles are spaced by 0.1 Gyr and the large circles by 1 Gyr. A range 
of redshifts are also indicated by the numbers along the track. The 
dashed line is the track traced out by the star-bursting galaxy M82. 
The diamonds indicate intervals of 1 Gyr, equivalent to the large 
circles for the SBc template. The solid line indicates the color-cut 
used to separate high- and low-redshift sources — sources above the 
line are mostly at z > 1.2 while sources below t he line are mostly at 
z < 1.2. The choice of this line is described in|Dcvlin ct alj J2009I . 
supplement) . We use Hp = 70.5k ms- 1 Mpc" 1 , S!m = 0.274, and 
£7a = 0.726 (Hinshaw et al. 2009) to tabulate the lookback times. 

millimeter galaxies (jPope et al .1 12008). is also shown. Its 
IRAC colours as a function of redshift are very similar 
to the SBc galaxy. 
We adopt the line 

([3.6] - [4.5]) = 

0.0682 x ([5.8] - [8.0]) - 0.075, (8) 

where the quantities in square brackets are AB magni- 
tudes in the IRAC bands, to divide the sources into two 
redshift bins. We have checked this cut for the 4242 
sources in the MUSIC catalog for which spectroscopic or 
photometric redshifts are available, finding that sources 
lie at redshifts higher and lower than z = 1.2 above and 
below the line, respectively. The cut i s remarkably sharp , 
with only 15% cross-contamination (Devlin et al.l [2009L 
supplem ent). This color- color cut is similar to the cut 
used by lYun et all (|2008l Figure 2) to identify submil- 
limeter galaxy counterparts. 

The differential contributions to the CIB from these 
"high" and "low" redshift sources as a function of 24 fim 
flux density are shown in Figure [H We also show the 
corresponding curve for the total CIB (using the entire 
FIDEL catalog) for reference. Specifically, the stacks 
have been split into logarithmic bins in flux density and 
the total stacked intensity is determined for each bin, 
divided by the linear bin widths. We therefore plot 
(Al I AS24) 5*24 as a function of S24, where Al is the to- 
tal contribution to the CIB at a given BLAST frequency 
from the sources in each bin and AS24 is the linear bin 
width. 

The bulk of the CIB at 250 and 350 fim is clearly pro- 
duced by sources with S24 > 60 /^Jy, while the low-flux 
density drop-off is less clear at 500 fim. This drop-off 
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Fig. 8. — The completeness-corrected differential distribution of 
contribution to the CIB versus flux density of the corresponding 
24 /jm galaxies. The result for the uncorrected differential stack 
is also shown for reference with a grey solid line to emphasize the 
fact that completeness only affects the two lowest bins. The 24 (im 
s ourc es are divided into high and low redshift bins, as described in 
§ 14.21 The binning in 24 (im flux density is the same as in Figure [5] 
with the 3 highest-flux density bins suppressed due to the small 
number of sources. The lowest 24 fim flux density bin is comprised 
of tentative detectio ns, a nd so is possibly less certain than the error 
bars indicate (see § 14.51 1. The low- and high-redshift curves have 
been displaced horizontally by —3 and +3%, respectively, for visual 
clarity. 

at low 24 /im flux densi ty is due in large p art to the 
number counts at 24 /im (jPapovich et al.ll2004l Figure 2), 
which are steepest at ~ 200-1000 /iJy and roll over sig- 
nificantly at lower flux densities. There appears to be 
CIB intensity missing from our measurements due to the 
flux density limit of the 24 /im catalog, in the sense that 
one would not expect strictly zero BLAST intensity in 
the next lower bin missing from each panel in Figure [H 
We expect that this missing component is small. The 
24 /im number counts below ~ 20 /iJy are uncertain, but 
since the 2 a upper limits derived from FIRAS are not 
much above the values derived in this analysis (Table [1]), 
we conclude that there is no hidden population of very 
faint 24 /jm sources that are important to the CIB at the 
BLAST wavelengths. 

Two trends in redshift are clear: first, that the fraction 
of the total CIB due to high redshift galaxies increases 
with wavelength through the BLAST bands (compare 
the relative heights of the dot-dashed and dashed curves 
in each panel); second, in all BLAST bands, the rela- 
tive contribution to the background produced by high- 
redshift sources (i.e. observed to be colder) increases to- 
ward fainter 24 /im flux densities (compare the heights of 
the dot-dashed and dashed curves as a function of S24). 
This analysis is consistent with our interpretation of Fig- 
ure [6] described in § 14.11 

Figure [9] shows the fraction of the CIB produced by the 
high-redshift sample as a function of wavelength. The 
curves are predicti ons from the phenom enological evo- 
lutionary model of Valiantc ct al. (2009). The different 
curves show the ratios of total intensities from galax- 
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Fig. 9. — The fraction of the 24 /^m-driven CIB which comes 
from high redshift galaxies, as defined by the IRAC color-color 
cut shown in Figure[7] is plotted against wavelength for the three 
BLAST bands and SCUBA. Over-plotte d are prediction s from the 
phenomenological evolutionary model of Valiantc ct al. (2009), for 
a range of redshift cuts. The solid, dashed, and dash-dotted curves 
are for a 24 /im catalog complete above 524 = 30 /^Jy. The dotted 
curve shows how to the z c = 1.4 curve changes for a S24, > 10 fijj 
catalog. z c = 1.2 is clearly not a good fit to the data. We hy- 
pothesize that our estimate of z c corresponding to the redshift cut 
of § 14.21 is biased low due selection effects. Note that the BLAST 
uncertainties are highly correlated (Table [TJ. 

ies brighter than S24 = 30 /iJy (solid, dashed, and dot- 
dashed curves) and S24 = 10/iJy (dotted curve) from 
galaxies located at z > labeled. We note that 

the curve for z c = 1.4 is a substantially better fit than 
the curve for z c — 1.2, despite the fact that photometric 
and spectroscopic redshifts in the MUSIC catalog sug- 
gest that the IRAC color-color cut we have adopted cor- 
responds to z c = 1.2. Perhaps the subset of sources with 
redshifts from this training set is biased, and the galaxies 
dominating the stacks do, in fact, lie at slightly higher 
redshifts than the models suggest. W e note that a sim - 
ilar IRAC color-color cut proposed bv lYun et al.l (2008) 
is designed to identify counterparts to 850 /im-selected 
galaxies that lie at redshifts predominantly z > 2. Al- 
ternatively, there are parameters governi ng the density 
and lu minosity evolution of galaxies in the lValiante et all 
(2009) model that could mimic this effect. 

4.3. Galaxies hosting an AGN 

We attempt to select for active galactic nuclei (AGN) 
in the FIDEL sample based on color-color cuts. We 
attempted to use th e power-law selection suggested by 
iDonlev et alj (|2008f ). but found that the majority of 
sources, even those that look like a power law through 
the IRAC bands, are rejected by the goodness-of-fit cri- 
terion, due to the small relative errors in the IRAC flux 
densities. The majority of sources that met both the 
power law index and goodness-of-fit criteria were those 
with very poorly determined 5.8 and 8.0 /xm flux densi- 
ties, and thus landed all over the color-color plane. In- 
ste ad, we make a cut o n colors that encompasses all of 
the IDonlev et al.l (|2008f) sources with minimal contami- 
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Fig. 10. — The average BLAST flux density of 24 (im sources 
classified as AGN (solid lines) and for the full FIDEL list (dashed 
line). The classifi cation i s a col or-color cut based on the power law 
seletion of Donlev ct al. (2003). The lowest flux density bins are 
based on tentative 24 fiva detections, and are therefore biased low 
(see § 14.51 1. The curves are shifted horizontally by +5 and —5% for 
clarity. 

nation from non-AGN: 

y> 1.21 X- 0.30 (9) 
y< 1.21 x + 0.22 (10) 
y > -0.83 x + 0.16 (11) 

where y = log(S 8 ,o/S 4 . 5 ) and x = log (S 5 . 6 /S 3 . 6 ). We 
find that 480 FIDEL sources meet this criteria, which 
is similar to the expected number of X-ra y-detected 
AGN , based on the source density reported by (|Luo et al.l 
2008). The average BLAST flux density at the positions 
of the IRAC color-selected sources is 7.5 ±1.1, 7.4 ± 0.9, 
and 6.0 ± 0.8 mJy at 250, 350, and 500 /im, respectively. 
These values are significantly larger than for the full FI- 
DEL list, 4.6 ± 0.2, 3.6 ± 0.2, and 2.3 ± 0.2 mJy. 

Figure [TO] shows the average BLAST flux density 
binned by 24 fj,m flux density for the FIDEL sources la- 
beled as AGN (solid line) and for the full set (dashed 
line). We see that across all BLAST bands, the AGN 
are brighter than average at all but the greatest 24 /im 
flux densities, with the ratio (averaged over all 24 /im 
flux densities) increasing from 1.6 to 2.6, from 250 to 
500 /im. The fact that this fraction increases at longer 
wavelengths (under the assumption that higher-rcdshift 
sources should appear cooler) is puzzling since the space 
densities of both AGN and the ultra-luminous star- 
forming galaxies select ed at 850 /im grow quickly to red- 
shifts z ~ 2-2.5 (e.g. iChapman et alJ [20051 : iWall et~aTl 
2005, 2008). Either the dust in star-forming galaxies is 
generally warmer than in the hosts of AGN, or the to- 
tal background produced by purely star-forming galaxies 
occurs at lower redshifts. This latter explanation is more 
likely, since it is now believed that the bulk of the total 



stellar mass in the Universe was formed in more numer- 
ous, less luminous galaxies at lower redshifts than those 
detected in longer- waveleng th submillimeter surveys (e.g. 
Perez- Gonzalez ct al. 2008). 

For comparison, we also stack on catalogs of AGN se- 
lected by the Chandra X-ray Obs ervatory and by o ptical 
variability. The X-ray sample of ILuo et alj (|2008| ) con- 
tains 462 sources, 164 of which (ignoring upper limits) 
fall in the region of the i?-band vs. soft X-ray flux plot 
that is thought to select for AGN (|Luo et al.l '2008. Fig- 
ure 12). The average BLAST flux density at the posi- 
tions of the X-ray-selected AGN is 7.3±1.9, 7.5±1.5, and 
5.3±1.3mJy at 250, 350, and 500 /im, respectively, simi- 
lar to the values listed above. We note that the stack 
on the full list of 462 is marginally larger at 250 /jm 
(9.4 ± 1.2 mJy), although nearly identical at 350 and 
500 /xm (7.6± 1.9 and 5.1±0 8mJy ). The optical variabil- 
ity sample of iTrevese et al.l (|2008D contains 132 objects, 
with average flux densities of 9.1 ± 2.1, 7.1 ± 1.7, and 
4.4 ± 1.4mJy at 250, 350, and 500 /im, respectively. 

Our measurements are tantalizing since it is one of the 
major goals in modern cosmology to establish an evo- 
lutionary connection between the formation of the most 
massive galaxies and the growth of black holes. How- 
ever, we warn the reader that the color-color selection 
described above may be significantly contaminated by 
star-forming galaxies, which would introduce strong bi- 
ases such that the observed trends here may be unrelated 
to the underlying evolution in the AGN and star-forming 
galaxy populations. A catalog of AGN which is complete 
and free of interlopers would allow us to improve our un- 
derstanding of the relative importance of AGN in the 
formation of the CIB. 

4.4. BzK Galaxies 

Dad di et al.l (|2004f ) have shown that selecting galaxies 
that are dimmer in z-band than they are in B and K s 
yields a list which contains many star forming galaxies. 
Specifically, BzK galaxies are those for which 

A BzK = (z - K) AB - {B - z) AB > -0.2, (12) 

which tends to select star forming galaxies at z > 
1.4. Their associa tion with ULIRGs has bee n studied 
(jDaddi et al.ll2005D . and lHartlev et all (|2008l ) conclude 
from their mild clustering properties that they are asso- 
ciated with dark matter halos. 

We use photometry from the 1548 FIDEL sources with 
MUSIC catalog identifications to calculate A Bz k- There 
are 388 BzK galaxies in the region, while the BLAST 
5 a catalog contains only a handful of entries, so we 
know that the BzK criterion is not selecting BLAST 
sources in general. Separating all 1548 FIDEL galax- 
ies into bins by A Bz k and stacking on the galaxies in 
each bin we find only a very weak dependence of mean 
BLAST flux density on A Bz k- BzK galaxies, i.e. sources 
with A Bz k > —0.2, are mildly brighter than others; the 
FIDEL/MUSIC BzK sources have average BLAST flux 
densities 1.3-1.7 times greater than for the full list, and 
produce 32 ±6, 34 ±7, and 42 ±11% of the total BLAST 
intensity at 250, 350, and 500 /im, respectively. It is clear 
that BzK samples make an important contribution to 
the CIB, but this is largely because they are so numer- 
ous. 
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For comparison, iTakagi et all (|2007l ) find that BzK 
galaxies produce 9% of the CIB at 850 /im, as measured 
by FIRAS. If we restrict our sample to galaxies with 
Ks < 22.9, the limiting magnitude of the ITakagi et alJ 
sample, we find that the BzK galaxies make up 12, 10, 
and 12% of the CIB at 250, 3 50, and 500 fim, respectively, 
comparable to Ta kagi et alJ 

4.5. Spurious Detections in FIDEL Catalog 

As described in § l2.2.1i detections in the FIDEL cat- 
alog below 30 /zJy are considered tentative. This affects 
the lowest flux density bin in Figures El and [101 Since 
Figures [5] and [10] are average flux densities, the effect 
of spurious sources (which are presumably uncorrelated 
with the BLAST maps) is to lower the estimate, so the 
lowest flux density bin is probably biased low. The un- 
corrected curves in Figure [8] are sums over the FIDEL 
sources, so spurious sources do not affect the measure- 
ment, and there is no bias. Spurious sources would af- 
fect the completeness correction, however, and thus the 
curves labelled "Full FIDEL" should also be regarded as 
biased low. In addition, since the quoted uncertainties 
take into account the number of sources included in the 
stack, the error bars are likely smaller than they should 
be. 

5. CONCLUSIONS 

For the first time, we trace the full submillimcter in- 
tensity of the far-IR background using a specific cata- 
log of 24 /im-selected galaxies. We use the properties of 
these 24 galaxies to analyze the composition of the 
CIB. We have determined that the average submillimeter 
flux density varies non-linearly with 24 /im flux density; 
these results require that the underlying averaged SED 
evolves with redshift. We make precise measurements 
of how the fraction of the CIB produced at higher red- 



shifts varies with observed wavelength through the sub- 
millimeter spectrum. These values are lower than models 
predict. We determine the total submillimeter intensity 
from unresolved AGN using an IRAC color-color selec- 
tion, finding weak evidence that they produce a non- 
negligible fraction of the CIB, and a relatively larger 
fraction in the longer-wavelength BLAST bands. This 
result could be demonstrating that the epoch of peak 
star-formation is at lower redshifts than the peak in AGN 
space density, which might be expected if most stars 
form in lower-luminosity, lower-redshift galaxies than the 
ultra-luminous star-forming galaxies selected at 850 /iin. 
However, this result depends strongly on the AGN selec- 
tion criterion which may contain many sources without 
an AGN, and which is also expected to be incomplete at 
higher redshifts. We also stacked at the positions of BzK 
galaxies, finding that they contribute 32-42% of the CIB 
in the B LAST bands, consi stent with the contribution 
found bv ITakagi et~aK (|200l at 850 /im when the depths 
of the optical catalogs are considered. 
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